import csv, sys
import ConfigParser
from ..sequence import Protein, Peptide
from ..io import readEntryFromNCBI
from ..misc import poissonParameters, poisson
from .. import Elements, Modifications, Aminoacids
from math import log
from scipy.integrate import inf
from re import finditer
from scipy.stats import poisson as ssPoisson
from zlib import crc32

_multiSpec = 1
_bothXMod = 2
_multiXLink = 4
_multiSeq = 8
_headTail = 16


def guessXonMatch(match, protTuple, spec, Tolerance, TolUnit):
	thisSpectrum = spec.descriptionSearch(match.source().description()).spectra()[0]

	s = 0
	e = 0
	posToTest = []
	guessed = 0
	aaList = Aminoacids.sections()
	
	if match.NPeptide().hasXMod():
		# we can set xstart as the position of xlinker
		given = match.NPeptide()
		toGuess = match.CPeptide()
		protein = protTuple[1]
		for (pos, mod) in match.NPeptide().modifications().items():
			if Modifications.getint(mod, 'xlinker') == 1:
				connections = Modifications.get(mod, 'links')
				s = pos
	else:
		# we can set xend as the position of xlinker
		given = match.CPeptide()
		toGuess = match.NPeptide()
		protein = protTuple[0]
		for (pos, mod) in match.CPeptide().modifications().items():
			if Modifications.getint(mod, 'xlinker') == 1:
				connections = Modifications.get(mod, 'links')
				e = pos
		
	for px in connections.split(','):
		if px in aaList:
			# by now we exclude the last aa... at time of writing the only xlinker
			# used prevents cleavage, so it is impossible to have the last aa modified
			# in a cleaved peptide... it can only be if the peptide is cterminal...
			if toGuess.isCTerminal(protein):
				seq = toGuess.sequence()
			else:
				seq = toGuess.sequence()[:-1]
			for pos in finditer(px, seq):
				posToTest.append(pos.start() + 1)
		elif px == 'np' and toGuess.isNTerminal(protein):
			posToTest.append(1)
		elif px == 'cp' and toGuess.isCterminal(protein):
			posToTest.append(toGuess.length())	

	if not posToTest:
		# there are no possible sites on other peptide
		return False
	elif len(posToTest) == 1:
		# there is only one site, we guess it is the right one
		guessed = posToTest[0]
	else:
		# there is more than one site
		# this is tough, because we now have to choose
		# we may decide this according to the highest 
		# number of ion found...
		# in posToTest I have index n... for each one I should check how many
		# bMass(n-1) and yMass(len-n) ions I have
		matchedIons = 0
		for i in posToTest:
			bIons = map(lambda x: toGuess.bMass(x), range(1, i - 1))
			yIons = map(lambda x: toGuess.yMass(x), range(1, toGuess.length() - i))
			# ok, I'm going weird with map and lambda and so on:
			try:
				nb = len(reduce(lambda x, y: x + y, map(lambda x: thisSpectrum.closeTo(x, Tolerance, TolUnit).peakMasses(), bIons)))
			except TypeError:
				nb = 0
				
			try:
				ny = len(reduce(lambda x, y: x + y, map(lambda x: thisSpectrum.closeTo(x, Tolerance, TolUnit).peakMasses(), yIons)))
			except TypeError:
				ny = 0

			if ny + nb < matchedIons:
				break
			else:
				matchedIons = ny + nb
		### BTW, if there are no matches... be strict, discard this!
		if matchedIons == 0:
			return False
		guessed = i

	if not s:
		s = guessed
	else:
		e = guessed

	match.setXStart(s)
	match.setXEnd(e)

	return True


class OMSSAMatch():
	""" Class containing any single omssa results match
	defined by its parameters
	"""
	def __init__(self):
		self._spectrumNumber = 0
		self._description = ''
		self._peptide = Peptide()
		self._mass = 0.0
		self._theoMass = 0.0
		self._charge = 1
		self._eValue = 0.0
		self._pValue = 0.0
		self._span = [0, 0]
		self._gi = 0
		self._proteinAccession = ''
		self._proteinDescription = ''
		self._crc32 = 0
	def setSpectrumNumber(self, value):
		""" set spectrum number
		"""
		self._spectrumNumber = int(value)
	def setDescription(self, value):
		""" set spectrum description
		"""
		self._description = value
	def setPeptide(self, value):
		""" set peptide string as result of OMSSA
		"""
		self._peptide.setSequence(value)
	def setPrecursorMass(self, value):
		""" set spectrum m/z
		"""
		self._mass = float(value)
	def setPrecursorTheoMass(self, value):
		""" set precursor Mr for this spectrum
		"""
		self._theoMass = float(value)
	def setSpan(self, value):
		""" Sets the span of the petided defined
		as start,end couple
		"""
		self._span = map(lambda x: int(x), value)
	def setGI(self, value):
		""" set identified protein gi
		"""
		self._gi = int(value)
	def setProteinAccession(self, value):
		""" set identified protein accession
		"""
		self._proteinAccession = value
	def setProteinDescription(self, value):
		""" set identified protein description
		"""
		self._proteinDescription = value
	def setCharge(self, value):
		""" set precursor charge
		"""
		self._charge = int(value)
	def setEValue(self, value):
		""" set e-value for this match
		"""
		self._eValue = float(value)
	def setPValue(self, value):
		""" set p-value for this match
		"""
		self._pValue = float(value)
	def setModificationsFromString(self, value):
		""" Sets the modifications for this match
		as a dictionary:
		mods{position:modification name,...}
		"""
		t = {}
		if value:
			for n in value.split(','):
				m = n.split(':')
				t[int(m[1])] = m[0]
		self._peptide.setModifications(t)
	def setAllWithCsvString(self, list):
		""" Sets any parameter for the match
		starting from a single csv results line 
		"""
		self.setSpectrumNumber(list[0])
		self.setDescription(list[1])
		self.setPeptide(list[2])
		self.setEValue(list[3])
		self.setPrecursorMass(list[4])
		self.setGI(list[5])
		self.setProteinAccession(list[6])
		self.setSpan([list[7], list[8]])
		self.setProteinDescription(list[9])
		self.setModificationsFromString(list[10])
		self.setCharge(list[11])
		self.setPrecursorTheoMass(list[12])
		self.setPValue(list[13])
		self._crc32 = crc32(list[0] + list[2] + list[10] + list[11])
	def crc32(self):
		return self._crc32
	def spectrumNumber(self):
		""" returns spectrum number
		"""
		return self._spectrumNumber
	def description(self):
		""" returns spectrum description
		"""
		return self._description
	def peptide(self):
		""" returns identified peptide
		"""
		return self._peptide
	def sequence(self):
		""" returns idetified peptide sequnece
		"""
		return self._peptide.sequence()
	def precursorMass(self):
		""" returns m/z
		"""
		return self._mass
	def precursorTheoMass(self):
		""" returns theoretical mass
		"""
		return self._theoMass
	def span(self):
		""" returns peptide span on identified protein
		"""
		return self._span
	def gi(self):
		""" returns protein gi
		"""
		return self._gi
	def proteinAccession(self):
		""" returns protein accession number
		"""
		return self._proteinAccession
	def proteinDescription(self):
		""" returns protein description
		"""
		return self._proteinDescription
	def charge(self):
		""" returns precursor charge
		"""
		return self._charge
	def eValue(self):
		""" returns match e-value
		"""
		return self._eValue
	def pValue(self):
		""" returns match p-value
		"""
		return self._pValue
	def modifications(self):
		""" returns peptide modifications as a dictionary
		"""
		return self._peptide.modifications()
	def length(self):
		return len(self.peptide())
	def error(self, unit = 'Da'):
		""" return error between obs and calculated
		precursor Masses, either in Dalton or in ppm
		"""
		if unit == 'Da':
			return self.precursorMass() - self.precursorTheoMass()
		elif unit == 'ppm':
			return (self.precursorMass() - self.precursorTheoMass()) / self.precursorTheoMass() * 1000000
	def precursorMz(self):
		""" return obs precursor m/z
		"""
		return (self.precursorMass() + (self.charge() * Elements.getfloat("H", "mmass"))) / self.charge()
	def bMass(self, ionIndex = 2):
		""" return the mass of b-ion with given ionIndex
		"""
		return self._peptide.bMass(ionIndex)
	def yMass(self, ionIndex = 1):
		""" return the mass of y-ion with given ionIndex
		"""
		return self._peptide.yMass(ionIndex)
	def aMass(self, ionIndex = 2):
		""" return the mass of a-ion with given ionIndex
		"""
		return self.bMass(ionIndex) - 27.9949

		
class OMSSAResults():
	""" Class that defines a single OMSSA results file
	made of a list of OMSSMatch
	"""
	def __init__(self):
		self._matchList = []
		self._proteins = {}
		self._crclist = set()
	def parseOMSSACsv(self, csvfile):
		""" Parse the csvfile and fills the class properties
		"""
		d = csv.excel
		d.lineterminator = '\n'
				
		try:
			fh = open(csvfile, 'r')
		except IOError:
			sys.stderr.write('Missing file ' + csvfile + '\n')
			sys.exit(1)
		for x in csv.reader(fh, dialect = d):
			if x[0] != 'Spectrum number':
				n = OMSSAMatch()
				n.setAllWithCsvString(x)
				if n._crc32 in self._crclist:
					continue
				self._crclist.add(n._crc32)
				self._matchList.append(n)
		fh.close()
	def addMatch(self, match):
		""" add a match as OMSSAMatch object to the
		match list
		"""
		if isinstance(match, OMSSAMatch):
			if match.crc32() not in self._crclist:
				self._matchList.append(match)
				self._crclist.add(match.crc32())
		else:
			pass
	def removeMatchAtIndex(self, index):
		""" remove match at given index from list and returns it
		"""
		self._matchList.pop(index)
	def removeLastMatch(self):
		""" remove last match from list and returns it
		"""
		self._matchList.pop()
	def removeFirstMatch(self):
		""" remove first match from list and returns it
		"""
		self.removeMatchAtIndex(0)
	def removeMatch(self, match):
		""" remove given match from list
		"""
		try:
			self._matchList.remove(match)
		except ValueError:
			None
	def items(self):
		""" return an iterator over matches
		"""
		return iter(self._matchList)
	def matches(self):
		""" return the complete list of matches
		"""
		return self._matchList
	def belowEValue(self, threshold):
		""" Returns a new OMSSAResults instance with
		all matches below a given e-value
		"""
		threshold = float(threshold)
		n = OMSSAResults()
		n._matchList = filter(lambda x: x.eValue() <= threshold, self._matchList)
		return n
	def belowPValue(self, threshold):
		""" Returns a new OMSSAResults instance with
		all matches below a given p-value
		"""
		threshold = float(threshold)
		n = OMSSAResults()
		n._matchList = filter(lambda x: x.pValue() <= threshold, self._matchList)
		return n
	def length(self):
		""" Returns the number of matches for this results file
		"""
		return len(self._matchList)
	def precursorMassSearch(self, mass, tolerance = 0.5, unit = 'Da'):
		""" Return a new OMSSAResults object with all the spectra
		with a precursor Mr equal to the mass provided, within
		a certain tolerace
		"""
		t = OMSSAResults()
		mass = float(mass)
		tolerance = float(tolerance)

		if unit == 'Da':
			None
		elif unit == 'ppm':
			tolerance = mass / 1000000 * tolerance
		else:
			return None
			
		t._matchList = filter(lambda x: abs(x.precursorMass() - mass) <= tolerance, self._matchList)		
		return t
	def precursorMzSearch(self, mass, tolerance = 0.5, unit = 'Da'):
		""" Return a new OMSSAResults object with all the spectra
		with a precursor m/z equal to the mass provided, within
		a certain tolerace
		"""
		t = OMSSAResults()
		mass = float(mass)
		tolerance = float(tolerance)

		if unit == 'Da':
			None
		elif unit == 'ppm':
			tolerance = mass / 1000000 * tolerance
		else:
			return None
			
		t._matchList = filter(lambda x: abs(x.precursorMz() - mass) <= tolerance, self._matchList)		
		return t
	def precursorTheoMassSearch(self, mass, tolerance = 0.5, unit = 'Da'):
		""" Return a new OMSSAResults object with all the spectra
		with a precursor calculated Mr equal to the mass provided, within
		a certain tolerace
		"""
		t = OMSSAResults()
		mass = float(mass)
		tolerance = float(tolerance)

		if unit == 'Da':
			None
		elif unit == 'ppm':
			tolerance = mass / 1000000 * tolerance
		else:
			return None

		t._matchList = filter(lambda x: abs(x.precursorTheoMass() - mass) <= tolerance, self._matchList)		
		return t
	def chargeSearch(self, charge):
		""" Return a new OMSSAResults object with all the spectra
		with a certain charge value
		"""
		t = OMSSAResults()
		charge = int(charge)
		t._matchList = filter(lambda x: x.charge() == charge, self._matchList)		
		return t
	def descriptionSearch(self, string):
		""" Return a new OMSSAResults object with all the spectra
		with given text into spectrum description. Useful for filtering
		selected spectra...
		"""
		t = OMSSAResults()
		t._matchList = filter(lambda x: string in x.description(), self._matchList)
		return t
	def checkPeptideConsistency(self, MinLength = 3):
		""" Returns a new OMSSAResults object in which only peptides
		that are consistent with xlink event are present.
		Also, peptides are filtered on length basis.
		"""
		filtered = OMSSAResults()
		XlinkerList = filter(lambda x: Modifications.getint(x, 'xlinker') == 1, Modifications.sections())

		for m in self.items():
			
			xcount = 0	# counter for x-modifications in results
			
			# defining NPeptide and CPeptide now is only useful to understand
			# if the x-mod is on the last aminoacid of single peptides
			
			NPeptide = map(lambda x: x.split(':'), m.proteinDescription().split())[0][0]
			if m.sequence().find(NPeptide) == 0:
				CPeptide = map(lambda x: x.split(':'), m.proteinDescription().split())[1][0]
			else:
				CPeptide = map(lambda x: x.split(':'), m.proteinDescription().split())[0][0]
				NPeptide = map(lambda x: x.split(':'), m.proteinDescription().split())[1][0]
		
			if len(NPeptide) < MinLength or len(CPeptide) <= MinLength:
				continue
		
			if NPeptide + CPeptide == m.sequence() or CPeptide + NPeptide == m.sequence():
				# x-peptide sequence should be a combination of single peptides
				# this is the first, slowest test...
				for (pos, mod) in m.modifications().items():
					if mod in XlinkerList:
						# first of all check that modification is not on
						# C-terminal of one of single peptides...
						# this check should be made atomic because I'm not sure
						# every xlinker needs this
						
						# pos, the position, is relative to x-peptide, we should map it
						# on single peptides.
						if pos == len(NPeptide) or pos == m.peptide().length():
							# len(NPeptide) corresponds to the length of left peptide
							# pos is given starting from 1, so that if the modification is on 
							# the C-terminal of Left peptide, pos == len(NPeptide)
							# for the same reason, pos == len(m.peptide()) means that modification
							# is on C-term of CPeptide...
							# in this case we set xcount to a value greater than 2 and get ot
							# this loop
							xcount = 2
							break # breaks this loop of modifications...
						
						# after all count how many xlinkers we have
						# they should not be more than one per peptide
						xcount += 1
				if xcount == 1:
					# we should find just one x-mod for each peptide
					filtered.addMatch(m)
		return filtered
	def proteins(self, database):
		""" Returns a dictionary of protein objects for proteins in this
		object instance
		"""

		if len(self._proteins) > 0:
			# proteins have been already calculated
			return self._proteins

		for m in self.items():
			# get accession for proteins
			# suppose we have a lot of items... we may first collect a limited number of accession
			# numbers and then instanciate proteins as dictionary values.
			# if we instanciate objects as we find accession numbers
			# we perform much more I/O
			for p in map(lambda x: x.split(':'), m.proteinDescription().split())[0][1].split(','):
				self._proteins[p] = None
			for p in  map(lambda x: x.split(':'), m.proteinDescription().split())[1][1].split(','):
				self._proteins[p] = None
		
		for p in self._proteins.keys():
			# as we said above, we prefer to iterate now that we have a limited number
			# of accessions
			self._proteins[p] = Protein(readEntryFromNCBI(database, p))
		return self._proteins
	def guessXLinks(self, spectra, Tolerance = 0.5, TolUnit = 'Da'):
		""" Guess the xlinks and returns a XLinkList object
		"""
		XList = XLinkList()
		for m in self.items():
		
			proteinList = [[],[]]
			for p in map(lambda x: x.split(':'), m.proteinDescription().split())[0][1].split(','):
				proteinList[0].append(self._proteins[p])
			for p in map(lambda x: x.split(':'), m.proteinDescription().split())[1][1].split(','):
				proteinList[1].append(self._proteins[p])
		
			# what remains to do:
			# 1- map every x-link event in a protein framework
			# 2- collect up data and do some summary
			
			tempMatch = XLink()
			tempMatch.initWithMatch(m)
			if not tempMatch.orderedPeptide():
				# it could be necessary to reverse protein list, to be consistent with "left" and "right"
				# definitions...
				proteinList.reverse()
			
			# we can define here x-link start and end
			# later we can transform them in protein-relative coordinates...
			
			
			
			for t in [(n, c) for n in proteinList[0] for c in proteinList[1]]:
				
				# t is a tuple with a couple of proteins
				# it is not atomic because a peptide can be found more than once in
				# the protein sequence. We should create a new atomic XLink instance
				# each time a peptide is within its protein
				# this means: check how many times tempMatch.NPeptide() is within t[0]
				# and how many times tempMatch.CPeptide() is within t[1]
		
				# before appending the new object we need to apply another filter
				# it happens that I have N-N modifications but peptides are not
				# N-terminal to the protein... this means that we should skip
				# these... actually could be better asking Alessio....
		
		###	remove wrong xlinks according to modification type and its position
		###	Update: remove these can be wrong, especially if xlink chemical modification
		###			is the same on terminal aa and other aa... Just comment it by now
				# NPeptide:
				invalid = False
				for mod in tempMatch.NPeptide().modifications().values():
					type = Modifications.get(mod, 'modtype')
					modString = Modifications.get(mod, 'residues')
					if type == 'np' and not tempMatch.NPeptide().isNTerminal() and tempMatch.NPeptide().sequence()[0] not in modString:
						invalid = True
					if type == 'cp' and not tempMatch.NPeptide().isCTerminal() and tempMatch.NPeptide().sequence()[-1] not in modString:
						invalid = True
				# CPeptide:
				for mod in tempMatch.CPeptide().modifications().values():
					type = Modifications.get(mod, 'modtype')
					modString = Modifications.get(mod, 'residues')
					if type == 'np' and not tempMatch.CPeptide().isNTerminal() and tempMatch.CPeptide().sequence()[0] not in modString:
						invalid = True
					if type == 'cp' and not tempMatch.CPeptide().isCTerminal() and tempMatch.CPeptide().sequence()[-1] not in modString:
						invalid = True
		
				if invalid:
					continue
		### end of remove

		
				PeptideMappings = [[],[]]
		
				
				offset = 0
		
				# the following two loops build PeptideMappings list
				# in PeptideMappings[0] we have all NPeptide positions on NProtein
				# in PeptideMappings[1] we have all CPeptide positions on CProtein
				# we have then to generate all the couples, that with all protein couples
				# give the entire set of single, atomic xlinks
				# one step remains, that is: we need to correctly spot x-link position!!!!
				# we may want to do this before all these loops, because it should be the same
				# for all NewXLinks... we should calculate it for tempMatch and then replicate
				# values in all NewXLinks...
				
				# Here on we try to guess the aminoacid involved in x-link. For each atom one
				# peptide is explained by its modification, the other one not. 
				# To do this we should check xlink capabilities in Xlink definition
				
				if not guessXonMatch(tempMatch, t, spectra, Tolerance, TolUnit):
					continue
				
				# Here starts correct position mapping on each protein within each single
				# atomic x-link
				while True:
					found = t[0].sequence()[offset:].find(tempMatch.NPeptide().sequence())
					if found == -1:
						break
					PeptidePosition = offset + found
					offset = PeptidePosition + len(tempMatch.NPeptide().sequence())
					PeptideMappings[0].append(PeptidePosition)
		
				offset = 0
				while True:
					found = t[1].sequence()[offset:].find(tempMatch.CPeptide().sequence())
					if found == -1:
						break
					PeptidePosition = offset + found
					offset = PeptidePosition + len(tempMatch.CPeptide().sequence())
					PeptideMappings[1].append(PeptidePosition)
					
				for pt in [(np, cp) for np in PeptideMappings[0] for cp in PeptideMappings[1]]:
					
					NewXLink = XLink()
					NewXLink.initWithMatch(m)
					NewXLink.setSourceProteins(NProtein = t[0], CProtein = t[1])
					NewXLink.setXStart(tempMatch.xStart())
					NewXLink.setXEnd(tempMatch.xEnd())
					NewXLink.setNMap(pt[0] + NewXLink.xStart())
					NewXLink.setCMap(pt[1] + NewXLink.xEnd())
					NewXLink.setSpectrum(spectra.descriptionSearch(m.description()).spectra()[0])
					XList.appendXLink(NewXLink)
		return XList
				
		
class XLink():
	""" XLink atomic entity
	"""
	def __init__(self):
		self._NPeptide = Peptide()
		self._CPeptide = Peptide()
		self._NProtein = Protein()
		self._CProtein = Protein()
		self._source = OMSSAMatch()
		self._start = 0
		self._end = 0
		self._NMap = 0
		self._CMap = 0
		self._hasXMod = False
		self._pValue = 1.0
		self._calculatedPValue = False
		self._spectrum = None
		self._qValue = 1.0
	def initWithMatch(self, match):
		""" initialize object from another OMSSAMatch instance
		trying to extend it
		"""
		if not isinstance(match, OMSSAMatch):
			raise AttributeError, "Not a OMSSAMatch in value"
		
		self._source = match
		
		NPeptide = map(lambda x: x.split(':'), match.proteinDescription().split())[0][0]
		if match.sequence().find(NPeptide) == 0:
			CPeptide = map(lambda x: x.split(':'), match.proteinDescription().split())[1][0]
		else:
			NPeptide = map(lambda x: x.split(':'), match.proteinDescription().split())[1][0]
			CPeptide = map(lambda x: x.split(':'), match.proteinDescription().split())[0][0]

		self._NPeptide.setSequence(NPeptide)
		self._CPeptide.setSequence(CPeptide)
		
		# handle modifications...
		for (pos, mod) in self._source.modifications().items():
			if pos <= self._NPeptide.length():
				# position is on NPeptide
				self._NPeptide.addModification(pos, mod)
				if Modifications.getint(mod, 'xlinker') == 1:
					self._NPeptide.setXMod(True)
			else:
				# position is on CPeptide...
				pos = pos - self._NPeptide.length()
				self._CPeptide.addModification(pos, mod)
				if Modifications.getint(mod, 'xlinker') == 1:
					self._CPeptide.setXMod(True)
		
	def setSourcePeptides(self, NPeptide = None, CPeptide = None):
		""" set source peptides. Peptides can be strings but
		Peptide() object is preferred
		"""
		if isinstance(NPeptide, str):
			self._NPeptide.setSequence(NPeptide)
		if isinstance(Cpeptide, str):
			self._Cpeptide.setSequence(CPeptide)
	def sourcePeptides(self):
		""" return a list with both peptides part of this atomic
		xlink occurrence
		"""
		return [self._NPeptide, self._Cpeptide]
	def NPeptide(self):
		""" return n-peptide, on which the xlink starts
		"""
		return self._NPeptide
	def CPeptide(self):
		""" return c-peptide, xlink end peptide
		"""
		return self._CPeptide
	def setSourceProteins(self, NProtein = None, CProtein = None):
		""" set source proteins. Proteins have to be Protein() objects
		"""
		if NProtein == None:
			pass
		elif isinstance(NProtein, Protein):
			self._NProtein = NProtein
		else:
			raise AttributeError, "NProtein is not a Protein() object"

		if CProtein == None:
			pass
		elif isinstance(CProtein, Protein):
			self._CProtein = CProtein
		else:
			raise AttributeError, "CProtein is not a Protein() object"
	def sourceProteins(self):
		""" returns source Proteins as list
		"""
		return [self._NProtein, self._CProtein]
	def NProtein(self):
		""" return protein the NPeptide belongs to
		"""
		return self._NProtein
	def CProtein(self):
		""" return protein the CPeptide belongs to
		"""
		return self._CProtein
	def setXStart(self, p):
		""" set x-link position on NPeptide
		"""
		if p < 0 or p > self._NPeptide.length():
			raise IndexError, "Position is outside peptide"
		self._start = int(p)
	def setXEnd(self, p):
		""" set x-link position on CPeptide
		"""
		if p < 0 or p > self._CPeptide.length():
			raise IndexError, "Position is outside peptide"
		self._end = int(p)
	def xStart(self):
		""" return xlink position on NPeptide
		"""
		return self._start
	def xEnd(self):
		""" return xlink position on CPeptide
		"""
		return self._end
	def setNMap(self, p):
		""" set NPeptide position on NProtein
		"""
		self._NMap = int(p)
	def setCMap(self, p):
		""" set CPeptide position on CProtein
		"""
		self._CMap = int(p)
	def NMap(self):
		""" return NPeptide position on NProtein
		"""
		return self._NMap
	def CMap(self):
		""" return CPeptide postion on CProtein
		"""
		return self._CMap
	
	def orderedPeptide(self):
		""" returns True if NPeptide is also the first peptide
		in match protein definition, False otherwise
		"""
		
		if self._NPeptide.sequence() == map(lambda x: x.split(':'), self._source.proteinDescription().split())[0][0]:
			return True
		return False
	def source(self):
		return self._source
	def setSpectrum(self, s):
		""" associate a MS/MS spectrum to this atom
		"""
		self._spectrum = s
	def spectrum(self):
		""" return the spectrum associated to this atom
		"""
		return self._spectrum
	def setPValue(self, p):
		""" associate a p-value for this atom
		"""
		self._pValue = p
	def pValue(self):
		""" returns the p-value for this atomic occurrence
		calculated with a cumulative poisson distribution
		"""
		if self._calculatedPValue:
			return self._pValue
		else:
			self.calcPValue()
			return self._pValue
	def calcPValue(self, Tolerance = 0.5, Unit = 'Da', AllIons = False):
		""" Calculate the p-value for this XLink
		"""
		xladder = set()
		xam = set()
		xbm = set()
		nMean = 0
		cMean = 0
		xMean = 0
		nPeaks = 0
		cPeaks = 0
		xPeaks = 0
		supposedToBeIn = 0

#		filteredSpectrum = self.spectrum().filterOnIntensity(0.01 * self.spectrum().mostIntense().peakIntensities()[0])
		filteredSpectrum = self.spectrum()
		for i in range(1, 1 + len(self.NPeptide().sequence())):
			# iterate over the N-Peptide
			revI = len(self.NPeptide().sequence()) - i
			if i < self.xStart():
				# if we are in the position at n-terminal of the supposed
				# xlink start calculate the expected value of the poisson distribution [tLambda] and the
				# number of peaks matched [Kappa]
				(Lambda, Kappa) = poissonParameters(filteredSpectrum, self.NPeptide().bMass(i), Tolerance, Unit)
				nMean += Lambda		
				nPeaks += Kappa		
				supposedToBeIn += 1
				# add to xam set the mass of the y-ion after the xlink position
				# these ions include also the xlinker itself... they will be used later on to
				# create a "complex ion ladder"
				xam.add(self.NPeptide().yMass(revI))
			else:
				# we are after the xlink position. We have to get y-ion values
				(Lambda, Kappa) = poissonParameters(filteredSpectrum, self.NPeptide().yMass(revI), Tolerance, Unit)
				nMean += Lambda
				nPeaks += Kappa				
				supposedToBeIn += 1
				# add to xam the mass of the b-ion before the xlink position
				xam.add(self.NPeptide().bMass(i))

		for i in range(1, 1 + len(self.CPeptide().sequence())):
			# iterate over the C-Peptide and do the same as above
			revI = len(self.CPeptide().sequence()) - i
			if i < self.xEnd():
				(Lambda, Kappa) = poissonParameters(filteredSpectrum, self.CPeptide().bMass(i), Tolerance, Unit)
				cMean += Lambda
				cPeaks += Kappa				
				xbm.add(self.CPeptide().yMass(revI))
				supposedToBeIn += 1
			else:
				(Lambda, Kappa) = poissonParameters(filteredSpectrum, self.CPeptide().yMass(revI), Tolerance, Unit)
				cMean += Lambda
				cPeaks += Kappa				
				xbm.add(self.CPeptide().bMass(i))
				supposedToBeIn += 1

		for i in xam:
			for o in xbm:
				# create a mass ladder that includes complex ions
				xladder.add(i + o)

		for m in xladder:
			# for this ladder caluclate the poisson values too
			(Lambda, Kappa) = poissonParameters(filteredSpectrum, m, Tolerance, Unit)
			xMean += Lambda
			xPeaks += Kappa				
		
		# now we can calculate the overall p-value as results of a poisson distribution caluclated
		# with expected values and number of matches
		if AllIons:
			self._pValue = poisson(nMean + cMean + xMean, nPeaks + cPeaks + xPeaks)
		else:
			self._pValue = poisson(nMean + cMean, nPeaks + cPeaks, supposedToBeIn)
#			if nPeaks + cPeaks == 0:
#				self._pValue = 0
#			self._pValue = poisson(self.NPeptide().length() + self.CPeptide().length() - 2, nPeaks + cPeaks)
		
#		self._pValue = poisson(nMean + cMean, nPeaks + cPeaks)
#		self._pValue = poisson(nMean, nPeaks) * poisson(cMean, cPeaks) * poisson(xMean, xPeaks)
		# set the flag "_calculatedPValue" to 1, so that we do not have to recalculate it...
		self._calculatedPValue = True
	def setQValue(self, qValue):
		self._qValue = qValue
	def qValue(self):
		return self._qValue
		
		
class Interaction():
	""" A class that describes an xlink at protein level
	"""
	def __init__(self):
		self._NProtein = Protein()
		self._CProtein = Protein()
		self._Peptides = set()
		self._linkers = set()
		self._start = 0
		self._end = 0
		self._rank = 0
		self._gScore = 0.0
		self._xLinks = []
		self._calculatedGScore = False
		self._pepOrder = []
		self._rawXLMod = set()
	def setNProtein(self, p):
		""" set the start protein for this event
		"""
		self._NProtein = p
	def setCProtein(self, p):
		""" set the end protein for this event
		"""
		self._CProtein = p
	def setStart(self, start):
		""" set the start of the xl
		"""
		self._start = int(start)
	def setEnd(self, end):
		""" set the end of the xl
		"""
		self._end = int(end)
	def setRank(self, r):
		""" set rank
		"""
		self._rank = r
	def NProtein(self):
		""" returns the start (n-) protein
		"""
		return self._NProtein
	def CProtein(self):
		""" returns the end (c-) protein
		"""
		return self._CProtein
	def rank(self):
		""" returns the rank
		"""
		return self._rank
	def xLinks(self):
		""" returns the list of atomic xlink
		occurences for this event
		"""
		return self._xLinks
	def length(self):
		return len(self._xLinks)
	def start(self):
		""" return the aa position on start protein
		where the xlinker is bound
		"""
		return self._start
	def end(self):
		""" return the aa position on end protein
		where the xlinker is bound
		"""
		return self._end
	def linkers(self):
		""" returns a set of xlinker found in this event
		a set is required as there can be many different 
		linkers (i.e. isotopic variants)
		"""
		return self._linkers
	def peptides(self):
		""" retursn a set of peptides that have been found
		in this event
		"""
		return self._Peptides
	def numberOfSpectra(self):
		tmp = {}
		return len(tmp.fromkeys(map(lambda x: x.source().spectrumNumber(), self._xLinks)).keys())
	def addXLink(self, m):
		""" add an atomic xlink to this event and automatically
		update event rank, linker and peptides
		"""
		if not isinstance(m, XLink):
			return False
		if m in self._xLinks:
			# this should check for duplicates, actually checks that we aren't addint the
			# very same object... the check is done on memory address... below we also
			# check that we are not adding a real duplicate... this happens when we merge several
			# result files bringing duplicates
			return False
		if len(self._xLinks) == 0:
			# the event is empy, fill it with this
			self._NProtein = m.NProtein()
			self._CProtein = m.CProtein()
			self._start = m.NMap()
			self._end = m.CMap()
			self._xLinks.append(m)
#			self._NPeptides.add(m.NPeptide().sequence())
#			self._CPeptides.add(m.CPeptide().sequence())
			self._Peptides.update([m.NPeptide().sequence(), m.CPeptide().sequence()])
			self._linkers.update(filter(lambda x: Modifications.getint(x, 'xlinker') , m.source().modifications().values()))
			self._pepOrder = map(lambda x: x.split(':')[0], m.source().proteinDescription().split())
			self._rawXLMod.update(set(filter(lambda x: Modifications.getint(x[1], 'xlinker') , m.source().modifications().items())))
		else:
			# before going on, check that this is not a duplicate
			# a duplicate has the same coords, the same mods, the same spectrum number and the same scan number
#			for existingXLink in self._xLinks:
#				if m.source().spectrumNumber() == existingXLink.source().spectrumNumber() \
#				and m.source().description() == existingXLink.source().description() \
#				and m.NMap() == existingXLink.NMap() \
#				and m.CMap() == existingXLink.CMap() \
#				and m.NPeptide().modifications().items() == existingXLink.NPeptide().modifications().items() \
#				and m.CPeptide().modifications().items() == existingXLink.CPeptide().modifications().items():
#					return False
		
		
		
			# assign a numerical score depending on how xlinks are similar, i.e
			# 0: xlinks are different
			# 1: xlinks are equal in position and proteins
			# 2: xlinks are equal in position, proteins but differ in peptide sequences
			# 3: xlinks are equal in position and proteins but have different directions: N-P = P-N
			# 4: xlinks are equal in positions and proteins but have different x-modifications.
			# each condition is tested and values are summed up...
			# still I have to tune values...
			
			if set([m.NProtein().accession(), m.CProtein().accession()]) == set([self._NProtein.accession(), self._CProtein.accession()]) \
			and set([m.NMap(), m.CMap()]) == set([self._start, self._end]):
				# in this case we are asserting that the current match is already in the collection
				# at least for mapping in protein structure...


				searchSet = m.NPeptide().modifications().items() + m.CPeptide().modifications().items()

				testSet = set(filter(lambda x: Modifications.getint(x[1], 'xlinker') , searchSet))
				
				# there can be the same spectrum that is matched several times during OMSSA search
				# it happens that, after validation, you get duplicates...
				if m.source().spectrumNumber() in map(lambda x:  x.source().spectrumNumber(), self._xLinks):
					# this spectrum has not been seen yet, we turn on bit 1 in mask 
					if self._rawXLMod.intersection(testSet) != testSet:
						# same spectrum says that xlinker can be in two x-positions
						self._rank |= _bothXMod
				else:
					# different spectra for this event
					self._rank |= _multiSpec

				# in the end we add to the collection
				self._rawXLMod.update(testSet)
				self._xLinks.append(m) # add the xLink to this set
			else:
				return False
	
				
			testSet = set(filter(lambda x: Modifications.getint(x, 'xlinker') , m.source().modifications().values()))

			if self._linkers.intersection(testSet) != testSet:
				# this is another ideal situation, two different xlinkers (isotopic variants ideally)
				# are found binding the same aminoacids... turn on the 3rd bit
				self._rank |= _multiXLink
				self._linkers.update(testSet)

			testSet = set([m.NPeptide().sequence(), m.CPeptide().sequence()]) # this is a set() of all sequences for the xlink event
			if self._Peptides.intersection(testSet) != testSet:
				# the current xlatom has different sequences
				# suppose the current event is made by TRLYPNKIR and ASDFGKLMNR,
				# the current xlatom (m) may have LMRTRLYPNKIR ad ASDFGKLMNR or other
				# combinations... we turn on the 4th bit of the rank
				# to me if a event has 4th and 5th bit on can be called a really TRUE xlink!
				self._rank |= _multiSeq
				self._Peptides.update(testSet)

			if self._xLinks[0].source().sequence().find(m.NPeptide().sequence()) != 0:
				# this is the best case to me... it means that OMSSA was able to see ions from both
				# peptides: we can find TRLYPNKIRASDFGKLMNR and ASDFGKLMNRTRLYPNKIR
				# we turn on the 5th bit
				self._rank |= _headTail
				
		
		return True
	def hasMultipleSpectra(self):
		if self._rank & _multiSpec == _multiSpec:
			return 1
		return 0
	def hasBothXModifications(self):
		if self._rank & _bothXMod == _bothXMod:
			return 1
		return 0
	def hasMultipleXLinkers(self):
		if self._rank & _multiXLink == _multiXLink:
			return 1
		return 0
	def hasMultipleSequences(self):
		if self._rank & _multiSeq == _multiSeq:
			return 1
		return 0
	def hasHeadTail(self):
		if self._rank & _headTail == _headTail:
			return 1
		return 0
	def eValue(self):
		""" returns the combined e-value as P(e) iterated over
		each Xlink atom in this event
		"""
		elist = map(lambda x: x.source().eValue(), self._xLinks)
		return reduce(lambda x, y: x * y, elist)
	def gScore(self, spectrumSampling = {}):
		""" returns a log-likelihood score as 2 *  sum(log(p(i)/H(i))) for i in 
		Xlink atoms. p(i) is the p-value calculated, H(i) is the null probability
		null is simply the frequency of equal atoms in this event
		"""
		
		if self._calculatedGScore:
			return self._gScore
		
		
		# first build frequencies
		sampling = {}
		k = log(len(self._xLinks))
		G = 0.0


		for m in self._xLinks:
			# this is crap... but it works
			# It's a key to store in a dict each xlatom count... to distinguish between xlatoms
			key = m.source().sequence() + str(m.source().modifications())
			try:
				sampling[key] += 1
			except KeyError:
				sampling[key] = 1

			# the same but with spectra...
			try:
				spectrumSampling[m.spectrum().description()] += 1
			except KeyError:
				spectrumSampling[m.spectrum().description()] = 1
		
		
		# then calculate G-score as 2 * sum(ln(p/H))
		# Need to study a way to do it in the loop above
		
		for m in self._xLinks:
			key = m.source().sequence() + str(m.source().modifications())

			try:
#				G = G + sampling[key] * log(m.pValue() / (float(sampling[key]) / k))
#				G = G + log(m.pValue() / k)
#				G = G + sampling[key] * log(m.pValue() / sampling[key])
#				G = G + sampling[key] * log(spectrumSampling[m.spectrum().description()] * m.pValue() / (float(sampling[key]) / k))
				G += sampling[key] * (log(spectrumSampling[m.spectrum().description()]) + k + log(m.pValue()) - log(sampling[key]))
			except OverflowError:
				self._gScore = -inf
			except ZeroDivisionError:
				return 'You are an asshole'
			else:
				self._gScore = 2 * G
		self._calculatedGScore = True
		return self._gScore
		
class XLinkList():
	""" A class with a list of xlinks 
	and interactions. The atoms are the single spectra that are xlinks
	the events are physically mapped xlinks, many atoms can be part of the
	same event
	"""
	def __init__(self):
		self._xlList = []	# The list of atomic xlinks
		self._eList = []	# The list of xlink events
	def appendXLink(self, xlink):
		""" add a xlink
		"""
		self._xlList.append(xlink)
	def addInteraction(self, interaction):
		""" add an event
		"""
		self._eList.append(interaction)
	def calculatePValues(self, Tolerance = 0.5, unit = 'Da', AllIons = False):
		""" iterate over all xlinks atoms and calculate
		their p-value, given a tolerance
		"""
		map(lambda x: x.calcPValue(Tolerance, unit, AllIons), self._xlList)
	def errorFilter(self, threshold):
		""" Returns a new XLinkList object filtered by ppm threshold on precursor masses
		"""
		t = XLinkList()
		t._xlList = filter(lambda x: abs((x.source().precursorMass() - x.source().precursorTheoMass()) / x.source().precursorMass() * 1000000) <= threshold , self._xlList)
		return t
	def pValueFilter(self, threshold):
		""" Returns a new XLinkList object filtered by p-value threshold on sigle atoms
		"""
		t = XLinkList()
		t._xlList = filter(lambda x: x.pValue() <= threshold, self._xlList)
		return t
	def length(self):
		""" Returns the number of xlink atoms
		"""
		return len(self._xlList)
	def items(self):
		""" Iterate over the xlink atoms
		"""
		return iter(self._xlList)
	def interactions(self):
		""" Iterate over the events
		"""
		return iter(self._eList)
	def eLength(self):
		""" Returns the number or events
		"""
		return len(self._eList)
	def guessInteractions(self):
		""" Tries to assign each xlink in xlList to an event
		"""
		spectrumSampling = {}
		for m in self.items():
		
			scan = 0
			for x in self._eList:
				if not x.addXLink(m):
					scan += 1
		
			if scan == len(self._eList):
		
				self._eList.append(Interaction())
				self._eList[-1].addXLink(m)
		
			# correct spectrumSampling, this will be used later to correct on 
			# G-score caclulation... The idea is that the more a spectrum is matched
			# the less is valuable...
			
			try:
				spectrumSampling[m.spectrum().description()] += 1
			except KeyError:
				spectrumSampling[m.spectrum().description()] = 1
		
		map(lambda x: x.gScore(spectrumSampling), self._eList)
	def gScoreFilter(self, threshold):
		""" Returns a new XLinkList object filtered by gScore threshold on events
		"""
		t = XLinkList()
		t._xlList = self._xlList[:]
		t._eList = filter(lambda x: x.gScore() <= threshold, self._eList)
		return t
	def rankFilter(self, value, mode='gt'):
		""" Returns a new XLinkList object filtered by rank bitmask on events
		"""
		t = XLinkList()
		t._xlList = self._xlList[:]
#		t._eList = filter(lambda x: x.rank() >= bitMask, self._eList)
		if mode == 'gt':
			# keep if rank is higher than given value
			t._eList = filter(lambda x: x.rank() >= value, self._eList)
		elif mode == 'bit':
			# keep if rank is & with given value
			t._eList = filter(lambda x: x.rank() & value == value, self._eList)
		return t
	def specLengthFilter(self, value):
		""" filter off all events with a number of spectra less than given value
		"""
		t = XLinkList()
		t._xlList = self._xlList[:]
		t._eList = filter(lambda x: x.numberOfSpectra() >= value, self._eList)
		return t
		
	def sort(self, key = 'gScore'):
		""" Sort in place the event list based on the key
		"""

		if key == 'gScore':
			self._eList.sort(key = lambda obj:obj.gScore())

		if key == 'rank':
			self._eList.sort(key = lambda obj:obj.rank(), reverse = True)

	def correctByFDR(self, qThreshold = 0.05):
		""" filters in place the xlList with a threshold for FDR
		One needs to re-guess the events after this
		"""
		self._xlList.sort(key = lambda obj:obj.pValue())
		n = self.length()
		qValues = []
		for x in range(n):
			q = self._xlList[x].pValue() * n / (x+1)
			self._xlList[x].setQValue(q)
		self._xlList = filter(lambda x: x.qValue() <= qThreshold, self._xlList)
	def fdr(self, alpha = 0.05):
		""" returns calculated FDR for this analysis at a given alpha
		"""
		self._xlList.sort(key = lambda obj:obj.pValue())
		n = self.length()
		qValues = []
		for x in range(n):
			if self._xlList[x].pValue() > alpha * (x + 1) / n:
				return self._xlList[x].pValue()
	def setFixedModifications(self, ModList = ['carbamidomethyl C']):
		""" Sets fixed modifications. They should be passed as a list of modification
		titles. Moreover, this method should only be called once Npeptide and Cpeptide have been
		identified"""
		for xAtom in self.items():
			nseq = xAtom.NPeptide().sequence()
			cseq = xAtom.CPeptide().sequence()
			for modification in ModList:
				modtype = Modifications.get(modification, 'modtype')
				modresidues = Modifications.get(modification, 'residues').split()
				if modtype == 'aa':
					for x in range(len(nseq)):
						if nseq[x] in modresidues:
							xAtom.NPeptide().addModification(x + 1, modification)
					for x in range(len(cseq)):
						if cseq[x] in modresidues:
							xAtom.CPeptide().addModification(x + 1, modification)
				elif modtype == 'np':
					if xAtom.NPeptide().isNTerminal() and nseq[0] in modresidues:
						xAtom.NPeptide().addModification(1, modification)
					if xAtom.CPeptide().isNTerminal() and cseq[0] in modresidues:
						xAtom.CPeptide().addModification(1, modification)
				elif modtype == 'cp':
					if xAtom.NPeptide().isCTerminal() and nseq[-1] in modresidues:
						xAtom.CPeptide().addModification(len(nseq), modification)
					if xAtom.CPeptide().isCTerminal() and cseq[-1] in modresidues:
						xAtom.CPeptide().addModification(len(cseq), modification)
	def removedups(self):
		""" removes duplicated entries... they shouldn't be here but they are
		"""
		
		mapSet = set([])
		newList = []
		for xlink in self._xlList:
			sn = xlink.source().spectrumNumber()
			sc = xlink.source().charge()
			ss = xlink.source().peptide().sequence()
			sm = xlink.source()._rawModificationString
			k = str(sn) + ':' + str(sc) + ':' + ss + sm
#			npmap = xlink.NProtein().accession() + ':' + str(xlink.NMap())
#			cpmap = xlink.CProtein().accession() + ':' + str(xlink.CMap())
			seen = False
			if k in mapSet:
				seen = True
			else:
				mapSet.add(k)
			if not seen:
				mapSet.add(k)
				newList.append(xlink)
		self._xlList = newList[:]
		del mapSet
		del newList

		